
#ifndef AFTERMATH_NUMERICS_DSSE_COMPLEX_HPP_INCLUDED
#define AFTERMATH_NUMERICS_DSSE_COMPLEX_HPP_INCLUDED

#include "complex.hpp"
#include "dsse_math.hpp"

#include <emmintrin.h>
#include <mmintrin.h>
#include <smmintrin.h>

namespace aftermath
{
    namespace numerics
    {
        template <>
        struct complex<double>
        {
            typedef double data_type;

            complex() { }
        
            complex(double real, double imaginary)
            {
                this->m_data = _mm_set_pd(imaginary, real); 
            }
        
            complex& operator =(const complex& other) 
            { 
                this->m_data = other.m_data; 
                return *this; 
            };
            
            /** Warning: this operation is slow! */
            double get_real() const { return this->m_data.m128d_f64[0]; }
            /** Warning: this operation is slow! */
            double get_imaginary() const { return this->m_data.m128d_f64[1]; }

            void set_real(double value) { this->m_data = _mm_blend_pd(_mm_set_sd(value), this->m_data, 0x2); }
            void set_imaginary(double value) { this->m_data = _mm_blend_pd(this->m_data, _mm_set1_pd(value), 0x2); }
            
            complex operator +(const complex& other) { return complex(this->m_data + other.m_data); }

            complex operator -(const complex& other) { return complex(this->m_data + other.m_data); }

            complex operator *(const complex& other)
            {
                // this: (a, b), other: (c, d)
                __m128d tmp0, tmp1;
                complex result;
                tmp0 = _mm_shuffle_pd(other.m_data, other.m_data, 0x0); // (c, c)
                tmp1 = _mm_shuffle_pd(this->m_data, this->m_data, 0x1); // (b, a)
                result.m_data = _mm_mul_pd(this->m_data, tmp0); // (ac, bc)
                tmp0 = _mm_shuffle_pd(other.m_data, other.m_data, 0x3); // (d, d)
                tmp1 = _mm_mul_pd(tmp0, tmp1); // (bd, ad)
                result.m_data = _mm_addsub_pd(result.m_data, tmp1); // (ac - bd, bc + ad)
                return result;
            }

            complex operator /(const complex& other)
            {
                // this: (a, b), other: (c, d)
                __m128d tmp0, tmp1;
                complex result;
                tmp0 = _mm_shuffle_pd(other.m_data, other.m_data, 0x3); // (d, d)
                tmp1 = _mm_shuffle_pd(this->m_data, this->m_data, 0x1); // (b, a)
                result.m_data = _mm_mul_pd(this->m_data, tmp0); // (ad, bd)
                tmp0 = _mm_shuffle_pd(other.m_data, other.m_data, 0x0); // (c, c)
                tmp1 = _mm_mul_pd(tmp0, tmp1); // (bc, ac)
                result.m_data = _mm_addsub_pd(tmp1, result.m_data); // (bc - ad, ac + bd)
                result.m_data = _mm_shuffle_pd(result.m_data, result.m_data, 0x1); // (ac + bd, bc - ad)
                tmp0 = _mm_dp_pd(other.m_data, other.m_data, 0x33); // (cc + dd, cc + dd)
                result.m_data = _mm_div_pd(result.m_data, tmp0); // (ac + bd, bc - ad) / |cc + dd|^2
                return result;
            }

            complex& operator +=(const complex& other) { this->m_data += other.m_data; return *this; }

            complex& operator -=(const complex& other) { this->m_data -= other.m_data; return *this; }

            complex& operator *=(const complex& other)
            {
                // this: (a, b), other: (c, d)
                __m128d tmp0, tmp1;
                tmp0 = _mm_shuffle_pd(other.m_data, other.m_data, 0x0); // (c, c)
                tmp1 = _mm_shuffle_pd(this->m_data, this->m_data, 0x1); // (b, a)
                this->m_data = _mm_mul_pd(this->m_data, tmp0); // (ac, bc)
                tmp0 = _mm_shuffle_pd(other.m_data, other.m_data, 0x3); // (d, d)
                tmp1 = _mm_mul_pd(tmp0, tmp1); // (bd, ad)
                this->m_data = _mm_addsub_pd(this->m_data, tmp1); // (ac - bd, bc + ad)
                return *this;
            }

            complex& operator /=(const complex& other)
            {
                // this: (a, b), other: (c, d)
                __m128d tmp0, tmp1;
                tmp0 = _mm_shuffle_pd(other.m_data, other.m_data, 0x3); // (d, d)
                tmp1 = _mm_shuffle_pd(this->m_data, this->m_data, 0x1); // (b, a)
                this->m_data = _mm_mul_pd(this->m_data, tmp0); // (ad, bd)
                tmp0 = _mm_shuffle_pd(other.m_data, other.m_data, 0x0); // (c, c)
                tmp1 = _mm_mul_pd(tmp0, tmp1); // (bc, ac)
                this->m_data = _mm_addsub_pd(tmp1, this->m_data); // (bc - ad, ac + bd)
                this->m_data = _mm_shuffle_pd(this->m_data, this->m_data, 0x1); // (ac + bd, bc - ad)
                tmp0 = _mm_dp_pd(other.m_data, other.m_data, 0x33); // (cc + dd, cc + dd)
                this->m_data = _mm_div_pd(this->m_data, tmp0); // (ac + bd, bc - ad) / |cc + dd|^2
                return *this;
            }

            complex operator +(const double& right) { return complex(this->m_data + _mm_set1_pd(right)); }
            complex operator -(const double& right) { return complex(this->m_data - _mm_set1_pd(right)); }
            complex operator *(const double& right) { return complex(this->m_data * _mm_set1_pd(right)); }
            complex operator /(const double& right) { return complex(this->m_data / _mm_set1_pd(right)); }

            complex& operator +=(const double& right) { this->m_data = this->m_data + _mm_set1_pd(right); return *this; }
            complex& operator -=(const double& right) { this->m_data = this->m_data - _mm_set1_pd(right); return *this; }
            complex& operator *=(const double& right) { this->m_data = this->m_data * _mm_set1_pd(right); return *this; }
            complex& operator /=(const double& right) { this->m_data = this->m_data / _mm_set1_pd(right); return *this; }

            friend complex operator +(const double& left, const complex& right);
            friend complex operator -(const double& left, const complex& right);
            friend complex operator *(const double& left, const complex& right);
            friend complex operator /(const double& left, const complex& right);

        private:
            __m128d m_data;
            
            complex(__m128d& data) { this->m_data = data; }
        };

        complex<double> operator +(const double& left, const complex<double>& right) { return complex<double>(_mm_set1_pd(left) + right.m_data); }
        complex<double> operator -(const double& left, const complex<double>& right) { return complex<double>(_mm_set1_pd(left) - right.m_data); }
        complex<double> operator *(const double& left, const complex<double>& right) { return complex<double>(_mm_set1_pd(left) * right.m_data); }
        complex<double> operator /(const double& left, const complex<double>& right) 
        {
            complex<double> result;
            result.m_data = _mm_set_pd(-left, left);
            result.m_data *= right.m_data;
            result.m_data /= _mm_dp_pd(right.m_data, right.m_data, 0x33);
            return result; 
        }
    }
}

#endif // AFTERMATH_NUMERICS_DSSE_COMPLEX_HPP_INCLUDED
